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Numerical Integration 
1- Trapezoidal Rule 
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Figure . Trapezoidal rule. Figure. Composite trapezoidal rule. 


The composite trapezoidal rule. 
m-1 


I= J l= {f(x + 2 f(2) + 2 f(6) +۰° +2 f1) + f(0] 3 
1=1 


Example 
Suppose we wished to integrate the function trabulated the table 


below for f ( ا‎ over the Interval from x=1.8 to x=3.4 using n=8 


b 3.4 
Am=If (x)= ٣ (e“)dx 


XxX 1.6 | 1.8 2 2.2 2.4 | 2.6 | 28 3 3.2 | 3.4 3.6 | 8 


Solution 


%$---Trapezoidal Rule-------—--—-———-—- 
clc 
a=1.8; 
b=3.4; 
h=0.2; 
n= (b-a) /h 
£=0; 
x=2; 
for i=1:n; 
%c=a+ (i-1/2) *h; 
$fÊ=f+ (c^2+1); 
f= (ftexp (x) ) 
x=x+h; 
end 
Am_approx=h/2* (exp (a) +2*f+exp (b) ) 
syms +t 
Am exact=int (exp (t) ,1.8,3.4) 
error=Am exact-Am approx 
E_t= (error/ (Am_ approx+error) ) *100 


E_a=( (Am_ approx-Am exact) /Am_ approx) *100 


1 
1 
1 1 4 x X1 2با‎ ۸ 
X=a Xa Xa=b a b 


: Figure, Composite Simpson's 1/3 rule. 
Figure .. Simpson's 1/3 rule. ِ ٠ 


The composite Simpson’s 1/7 rule 
0 2 E RB 
flax e 1= [f(x) + 4 f(x) + 2 fa) +4 F(a) + °°° + 2 F(Xn-a) + 4 F(X) F FDIS 


A 


Example 
Suppose we wished to integrate the function using Simpson’s 1/4 rule and 


Simpson’s 34 rule the table below for f (x)=e* over the interval from x=1.8 


to x=3.4 using n=8 


b 3.4 
Am=Iff e ا‎ (e“)dx 


XxX 1.6 | 1.8 2 2.2 2.4 | 2.6 | 28 3 3.2 | 3.4 3.6 | 38 


Solution 


کے کے کے کڈ کے کے کے کے ھک کد کد ےد دک کد == = === SIMPSONS 1/3 Elê‏ == 
clc‏ 

a=1.8; 

b=3.4; 

h=0.2; 

n= (b-a) /h 

£=0; 

m=0; 
for x=2: (hth) :3.2; 

f= (ftexp (x)) 7; 

end 


for x=2.2: (hth) :3; 
m= (m+exp (x) ) ; 
end 
Am_approx=h/3* (exp (a) +4*f+2*m+texp (b) ) 
syms +t 
Am exact=int (exp (t) ,1.8,3.4) 
error=Am exact-Am approx 


|د 


E_t= (error/ (Am_approxterror) ) *100 
E_a=( (Am_approx-Am exact) /Am_ approx) *100 


3-Simpson ’s 3/8 rule 


The composite Simpson’s 3/4 rule 
b 
f(x) dx # I1 = [f(xı) +3 flax) + 3f) + 2f) + ---- + 3 f(x) + 3 f(xn-1) + f(%n)] 


h=0.2; 


for x=2:h:2+h; 
f=f+exp (x) 


x=x+h; 


m=exp (x) ; 


for x=2.6:h:2.6+h; 
f=f+texp (x) ; 


end 


x=x+h; 
m=m+texp (x) ; 
x=x+h; 


f=f+texp (x) ; 


syms t 
Am exact=int (exp (t) ,1.8,3.4) 
error=Am exact-Am approx 
E_t= (error/ (Am_approx+terror) ) *100 
E_a=( (Am_approx-Am_ exact) /Am_ approx) *100 


for x=2:h:3.2; 
switch x 
case {2,2.2} 
f=f+exp (x) 
case {2.4} 
m=exp (x) ; 
case {2.6,2.8} 
f=f+texp (x) ; 
case {3} 
m=mtexp (x) ; 
otherwise 


f=f+exp (x) ; 


syms +t 
Am exact=int (exp (t) ,1.8,3.4) 
pretty (Am exact) 
error=Am exact-Am approx 
pretty (error) 
E_t= (error/ (Am_approx+error) ) *100 
pretty (E_t) 
E_a=( (Am_approx-Am exact) /Am_ approx) *100 
pretty (E_a) 


4-Lagrange Interpolating Polynomial Method 


Lagrange’s interpolation method uses the formula 


,_ (E-KX-B)..(K-F) , (X-Ko(E-X)..(E-K) 
EN mme Og TNE E E II 
(Xx 7K, J(X > &))...(Xg > XJ (K, KJK > X,)...(X; ¬ Kg) 
(£-XxK)o)J(X-X,)...(X—Xg-1) )(ء‎ 
ا‎ 
(X, ¬ Kg)(Xq = K2)...(Xq = Xa-1) 
EXAMPLE 
Given the data points 
use Lagranges method to determine y al xX = 1. 
Solution 
2 lu — xax — X3) _ O1 - 201 —3) _1 
"(a — xax x3) (0-20-33) 3 
2 (x — x(x — X3) _ (O1— O —3) | 1 
“  )xa - x()xa -×xs( )2-0()2-3( 
ا‎ (x — xJ (xX — x<2) __ (O1— O1 —-2) 1 
° ` (ts — xl(xs =x) (@—-0G—2) 8 
OT 28 
Y= yl + yata + Jag = 2 +11 — E 4 

Lagrange’ s interpolation method ---------‏ ق 

clc 

x=1; 

Ssyms x 

ةة ةةة ةةة ةةة ةةة ةةة ت ك 

x1=0; 

x2=2; 

x<3=3; 

> تتت Sass‏ ك 

y0=7; 

y1=11; 

y2=28; 

ڑا چ ج چ چاخ چ ےھ کے ج جر ےھ چاچ چ کےا چ ج جخ چ کے چ کے چ چا چ چ کے و چ چ کے چ چ چ چ چ ج چاچ چ چ چ چ چ کے ك 


(x1-×<2( * )×1-×3( ( 
(x×2-x<1( * )×2-×3( ( 
(x3-x1( * )×3-×2( ( 


10= ) )×-×2( * )×-×3( ( 
11= ) )×-×1( * )×-×3( ( 
12= ) )×-×1( * )×-×2( ( 


E CC 
کس کا‎ 


Example 2 
Construct the polynomial interpolating the data by using 


Lagrange polynomials 


10= ) )×-×2( * )×-×3( ( / ) )×1-×2( * )×1-×3( ( 
11= ) )×-×1( * )>×-×3( ( / ) )×2-×1( * )×2-×3( ( 
12= ) )×-×1( * )×-×2( ( / ) )×3-×1( * )×3-×2( ( 


y=y0*10+y1*11+y2*12; 
collect (y) 
8 ت‎ 


clc 

syms x 

p=0; 

s=]1 1/2 31; 
f£=[3 -10 2]; 
n=length (s) ; 
for i=1:n; 


1=1; 
for j=1:n; 
ة٤ ([=~ہة)‎ 
1=((x-s(j))/(s(i)-s(j)))*1; 
end 
end 


p=l1l.*f(i)+p; 
end 


p=collect (Pp) 


Example 2 
Construct the polynomial interpolating the data by using 


Lagrange polynomials 


clc 


x=input (' enter value of x:') 
pP=0; 

s=]1 1/2 3]; 

f=[3 -10 2]; 

n=length (s) ; 

for i=1:n; 


1=1; 

for j=1:n; 
if ([=~ذ)‎ 

1=((x-s(j))/(s(i)-s(j)))*1; 
end 
end 
pP=1.*f(i)+tp; 
end 


Pr 
fprintf ('\n p(%53.3f£)=%5.4f£',x,P) 


> ےد ف فت فت % 
syms x‏ 
pP=0;‏ 
for i=1:n;‏ 
;1=1 
for j=1:n;‏ 
ر ([=~ة) £ 
1=((x-s(j))/(s(i)-s(j)))*1;‏ 
end‏ 
end‏ 
pP=1.*f(i)+tp;‏ 
end‏ 
p=collect (Pp)‏ 
ج د چ ج ےھ خ ھ کے ےھ ج کے چ % 


p = -283/10 -53/5 *x^2 + 419/10 *x 
enter value of x:5 
x =5 
م‎ )5.000(--0 


Example 3 


Find the area by lagrange polynomial using 3 nodes 


F() 6.04964 13.464 29.964 


Solution 


F1=13.464; 


10= ) )×-×2( * )×-×3( ( / ) )×1-×2( * )×1-×3( ( 
AO0=int (10,1.8,3.4) 
11= ) )×-×1( * )×-×3( ( / ) )×2-×1( * )×2-×3( ( 
A1l=int(11,1.8,3.4) 
12= ) )×-×1( * )×-×2( ( / ) )×3-×1( * )×3-×2( ( 
A2=int(12,1.8,3.4) 


F=FO0*AO+F1*A1+F2*A2 
collect (F) 


syms x 
format long 
p=0; 
s=[]1.8 2.6 3.4]; 
f= [6.04964 13.464 29.964]; 
n=length (s) ; 
for i=1:n; 
1=1; 
for j=1:n; 
ے٤ (ز[=~ہد)‎ 
1=))x-s)j))/(s(i)-s(j)))*l1l; 
end 
end 
A=int(l,s(1),s(n)) 
P=A*f (i) +p; 
end 


5-Mid Point Rule 


Example 


Find the mid point approximation for 
b 2 

Am= f f (x)dx= [(x^*2+1)dx 
a -1 


using n=O 


Solution 


%$---Mid Point Rule-------——————-—- 
clc 
a=-1; 
b=2; 
n=6; 
h= (b-a) /n; 
£=0; 
for i=1:n; 
c=a+ (i-1/2) *h; 
f=f+ (c^2+1); 


6- Taylor series 


A function f(x) which possesses all derivatives up to order n at a point x = x, can be expanded 
in a Taylor series as 


1"(xo) 2 ۴“) 
(xx) +... + 


21 1 


f(x) = f(xo) + F(xo)(x — xo) + (x xo) 


1! 


If x, = 0. reduces to 


f0 
f(x) = f(O)+PF(O)x+ ا‎ ea و‎ 


Example 


Compute the first three terms of the Taylor series expansion for the function 


y = f(x) = tank 
ata = [/4. 


Solution: 


The Taylor series expansion about point a is given by 
f.) = f(a) + 1(4) (x - )ك + )۾‎ - a)” + ا‎ aE 


and since we are asked to compute the first three terms, we must find the first and second deriva- 
tives of f(x) = tanx. 


From math tables, E $0 f (x) = sec x . To find f'(x) we need to find the first 


dl‏ . 2 2 ا 
derivative of sec x, so we let z = sec x. Then, usin q,Secx = secx ° tanx , We get‏ 


a E O e 
— = Dsecx—secX = 2secX(secX:‘ tAanXK) = 2sec Xx - tanx 
dx dx 

Next, using the trigonometric identity 


"i 2 
sec XxX = tan x+1 


and by substitution , we get, 


dz ۓ‎ T(x) = 2(tanx + 1)tanx 
dx 


Now, at point a = r/4 we have: 


1 
1 
سم 
+ 
ا 
1 
ت 
rb‏ 
2 
I‏ 
3 
اظ 
1 
ټ 
نا 
E.‏ 
+ 
فس 
س 
فم 
1 
ل 


fr) TT) AE) 
f (a)= f| 1 = tanl 7 / = 1 f(a)= F7 | 


and by substitution into (6.125), 


7 : ۱ <2 
f(x) =1+2|x- ادا‎ 
A A A 


We can also obtain a Taylor series expansion with the MATLAB taylor(f,n,a) function where f 
is a symbolic expression, n produces the first n terms in the series, and a defines the Taylor 
approximation about point a. 


The following MATLAB script computes the first 8 terms of the Taylor series 
expansion of y = f(x) = tanx about a = r/4. 


Taylor Seis ===‏ ج 
clc‏ 

a=pi/4; 

sym x 

y=tan (x) ; 


z=taylor (y,8,a) ; 
pretty (Zz) 


Example 


Express the function 


in a Maclaurin’s series. 


Solution: 


A Maclaurin’s series has the form, that is, 


0 n) 
f(x) = £)0) +f (0)x + O +... + 


For this function, we have f(t) = e and thus f(0) = 1. Since all derivatives are e, then, 
f(0) = f'(0) = f"(0) = ... = 1 and therefore, 


2 3 
e t t 
(f) = lt aT aT 


MATLAB displays the same result. 


syms +t 
fn=taylor (exp (t) ) ; 
pretty (fn) 


ك 


clc 
clear 
Ssyms x 
%p2=taylor (cos (x) ,7,Pi/4) 
format long 
x=pi/3 
true=cos (pPi/3) 
p1=1/2*2^ (1/2) -1/2*2^ (1/2) * (x-1/4*pi) 
et_1=( (true-p1) /true) *100 
p2=1/2*2^ (1/2) -1/2*2^ (1/2) * (x-1/4*pi) -1/4*2^ (1/2) * (x-1/4*pi) ^2 
et_2=( (true-p2) /true) *100 
p3=1/2*2^ (1/2) -1/2*2^ (1/2) * (x-1/4*pi) -1/4*2^ (1/2) * (x- 
1/4*pi) ^2+1/12*2^ (1/2) * (x-1/4*pi) ^3 
et_3=( (true-p3) /true) *100 
p4=1/2*2^ (1/2) -1/2*2^ (1/2) * (x-1/4*pi) -1/4*2^ (1/2) * (x- 
1/4*pi) ^2+1/12*2^ (1/2) * (x-1/4*pi) ^3+1/48*2^ (1/2) * (x-1/4*pi) ^4 
et_4=( (true-p4) /true) *100 
p5=1/2*2^ (1/2) -1/2*2^ (1/2) * (x-1/4*pi) -1/4*2^ (1/2) * (x- 
1/4*pi) ^2+1/12*2^ (1/2) * (x-1/4*pi) ^3+1/48*2^ (1/2) * (x-1/4*pi) ^4- 
1/240*2^ (1/2) * (x-1/4*pi) ^5 
et_5= ( (true-p5) /true) *100 
p6=1/2*2^ (1/2) -1/2*2^ (1/2) * (x-1/4*pi) -1/4*2^ (1/2) * (x- 
1/4*pi) ^2+1/12*2^ (1/2) * (x-1/4*pi) ^3+1/48*2^ (1/2) * (x-1/4*pi) ^4- 
1/240*2^ (1/2) * (x-1/4*pi) ^5-1/1440*2^ (1/2) * (x-1/4*pi) ^6 
et_6=( (true-p6) /true) *100 


ك2 
0 


1.047197551196598 
true = 


0. 500000000000000 


0.521986658763282 
et 1 = 
-4.3973317526541 


P2 = 
0.497754491403425 
et 2 = 
0.449101719315004 
P3 = 
0.499869146930044 
et 3 = 
0.02614 
p4 = 


0.500007550810613 
et 4 = 
-0.00151016212553 


0. 500000304000373 
et _ 5Š = 
-6.080007448616696e-005 
p6 = 
0.49999998779865 
et 6 = 
2.440274993187329e-006 


^2; 


^4- 


^4- 


clear 


x=-pi/2:0.1:pi/2; 


p1=1/2*2.^*(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) ; 
pP2=1/2*2.^*(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x-1/4*pi) . 
p3=1/2*2.^*(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 


1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) . 


pP4=1/2*2.^ (1/2) -1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 
1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) .^3+1/48*2.^ (1/2) * (x-1/4*pi) . 
pP5=1/2*2.^ (1/2) -1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 
1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) .^3+1/48*2.^ (1/2) * (x-1/4*Pi) . 


1/240*2.^ (1/2) * (x-1/4*pi) .^5; 


p6=1/2*2.^*(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 
1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) .^3+1/48*2.^ (1/2) * (x-1/4*pi) . 
1/240*2.^ (1/2) * (x-1/4*pi) .^5-1/1440*2.^ (1/2) * (x-1/4*pi) .^6; 


,'P6', 


tr=cos (x) ; 

plot (x,p1,'y') ,pause (1) 
hold on 

plot (x,p2,'r') ,pause (1) 
hold on 

plot (x,p3,'g') ,pause (1) 
hold on 

plot (x,p4, 'k') ,pause (1) 
hold on 

plot (x,p5, 'm') ,pause (1) 
hold on 

plot (x,p6,'b') ,pause (1) 
hold on 

plot (x,tr, 'c') ,pause (1) 
hold on 
legend('p1','p2','p3' 


,'P4','p5' 


grid 


Figures - Figure 1 
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File Edit Wiew 


û € ll 2 


ك 


clc 

clear 

x=-pi/2:0.1:pi/2; 

p1=1/2*2.^(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) ; 
p2=1/2*2.^*(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x-1/4*pi) .^2; 
p3=1/2*2.^*(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 
1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) .^3; 
p4=1/2*2.^(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 
1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) .^3+1/48*2.^ (1/2) * (x-1/4*pi) .^4; 
p5=1/2*2.^(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 
1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) .^3+1/48*2.^ (1/2) * (x-1/4*pi) .^4- 
1/240*2.^ (1/2) * (x-1/4*pi) .^5; 
p6=1/2*2.^(1/2)-1/2*2.^ (1/2) * (x-1/4*pi) -1/4*2.^ (1/2) * (x- 
1/4*pi) .^2+1/12*2.^ (1/2) * (x-1/4*pi) .^3+1/48*2.^ (1/2) * (x-1/4*pi) .^4- 
1/240*2.^ (1/2) * (x-1/4*pi) .^5-1/1440*2.^ (1/2) * (x-1/4*pi) .^6; 
tr=cos (x) ; 

subplot (331) 

plot (x,p1, 'y') ,pause (4) 

title ('p1') 

grid on 

subplot (332) 

plot (x,p2,'r') ,pause (4) 

title ('p2') 

grid on 

subplot (333) 

plot (x,p3,'g') ,pause (4) 

title ('p3') 

grid on 

subplot (334) 

plot (x,p4, 'k') ,pause (4) 

title ('p4') 

grid on 

subplot (335) 

plot (x,p5, 'm') ,pause (4) 

title ('p5') 

grid on 

subplot (336) 

plot (x,p6,'b') ,pause (4) 

title ('p6') 

grid on 

subplot (338) 

plot(x,tr,'c') 

title ('Tr') 

grid on 


Figures - Figure 1 
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1# Hêk a a" =| 8| DE aE 


Example 
Find first and second derivatives for F(<)=x^2+2x<+2 


Solution 


To find first and second derivatives of Pn (x) -----‏ تتت سس 


clc 
a=[1 2 3]; 
syms x 
p=a (1) ; 
for i=1; 
p=a (i+1) +x*p; 
end 
disp ('First derivative') 
p2=p+x*diff (Pp) 
disp (' Second derivative') 


p22=diff (p2) 


ج س ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت 


First derivative 


2+2*× 


Second derivative 


p22 = 


Example 
P4(x)=3x^4-10x^3-48x^2-2x+12 at r=6 deflate the polynomial 
with Horners algorithm Find P3(x). 


Solution 


Horner alogorithm------—-——————————-—‏ کد 

clc 

a=[]3 -10 -48 -2 127; 

r=6; 

b (1) =a (1) ; 

p=0; 

n=length (a) ; 

for i=2:n; 
b(i)=a(i)+r.*b(i-1); 

end 

syms x 

for i=1:n; 
p=p+b (i) *x^ (4-i) ; 

end 

disp ('P3 (x)=') 


3*x<^*3+8*x<^2-2 


Numerical Differentiation 
1-Finite Difference Approximations 


The derivation of the finite difference approximations for the derivatives of f (x) 
are based on forward and backward Taylor series expansions of f (x) about xX, 


such as 
۴ 1f a hٌ r hh r n 4) 0 
fix+ BR) = f(x + hf (x) + (x) + (x) + 0 (x) +... (a 
e E ¥ 1 س‎ HF ra e 
f(x — Mm = f(x) — hf (x) + F0 31 F(x) + 1 F&I (b 
fx +28) . )( +2hf(%) + 2 ET 2 pr + 2 ۶“ ( E (e 
FE-PE e f" , N > 


We also record the sums and differences of the series: 


: 5 3 FHF f ( : 
flx+ mM + flix — B) = 2 f(x) + FF f"(x) + f6 as (e) 
f(x+ BM -— f(x — RM = 2hf' (x E e ٠ (1) 
4h 
fx +20) + f(x = 2M) = 2 f(x) + ANP f(x) + e 4... (g) 
: SE 
flx +2R) - f(x 2M) = 4hf (x) + — f) + ° (h) 


First Central Difference Approximations 


The solution of Eq. (f) for f(x) is 


C+) — f(x —F h* 
f= f(x + ا‎ E 


Keeping only the first term on the right-hand side, we have 


1 flx+ Rh = f(x — 
')() ا لدد‎ ٠ إ‎ ©) 
ES 77 4 (hn) 
which is called the first central difference approximation for f'(x). The term O(h*]) 
reminds us that the truncation error behaves as lı. 
From Eq. (e) we obtain 
ıa _ JX+M -2f(0 + flx—h) FP mg, 
د(‎ ٠ ال ال‎ ٠ ل‎ F۴ ٭]‎ . 
f) ا‎ 
Or 

_ flx+h -2 f(x) + f(x — P) 


Ol 
2 + O(n) 


f" 


Central difference approximations for other derivatives can be obtained from 
Eqs. (ah) in a similar manner. For example, eliminating f(x) from Eqs. (f) and (h) 
and solving for f(x) yield 
2 flx+2mM -2 f(x + N) +2 f(x — mM — f(x — 2) 
E 2h3 


The approximation 


flix +2mM = 4f(x +h) +6 F(X) =4 f(x — mM + f(x = 2M) 
n^ 


+ 0) ( 
/)( = + 0) 
First Noncentral Finite Difference Approximations 


These expressions are called forward and backward finite difference 
approximations. 


Noncentral finite differences can also be obtained from Eqs. (a)—h). Solving 
Eq. (a) for f'(x) we get 


1 flx+ mM — f(x) 1 1 hn i î (4 fy 
وا و2 و‎ 


Keeping only the first term on the right-hand side leads to the first forward difference 
approximation 


i_J FG 
f= د‎ 


+ OR) 


Note that the truncation error is now O(N), which is not as good as the O(l*) error in 
central difference approximations. 


We can derive the approximations for higher derivatives in the same manner. For 
example, Eqs. (a) and (c) yield 


fx + 2R) - 2 f(x +R) + f(x) 
h2 


j)» = + O() 


Second Noncentral Finite Difference Approximations 


Finite difference approximations of O(h) are not popular due to reasons that will be 
explained shortly. The common practice is to use expressions of Oh). To obtain 
noncentral difference formulas of this order, we have to retain more terms in the 
Taylor series. As an illustration, we will derive the expression for] f(x). We start with 
Eqs. (a) and (c), which are 


۲ e MH TCR 
jflx+ Rh) = f(x) + hf ET eT iS (x) + <۰. 


۲ 2 ا‎ 4n rr 2h (4) f 
f(x + 2M) = f(x) + 2hf' (x) + 2h f" (x) + > f" (0) + = f) + ° 


We eliminate f"(x) by multiplying the first equation by 4 and subtracting it from the 
second equation. The result is 


5 2h 
FEMA = SFO RPO LE FO e 


Therefore, 
o — fix + 2h) +4 f(x + RM) — 3 f(x) fh“ f(x) + 
2h 3 
ol 
f) = rTM rte + 0)1 ( 


Tlus Equation is called the second forward finite difference approximation. 


EXAMPLE 


Use forward difference approximations of oh to estimate the first 


% derivative of 


fx = -0.1.*x.^4-0.15.*x.^3-0.5.*x.^2-0.25.*× +1.2 


solution 


% Use forward difference approximations to estimate the first 


% derivative of fx=-0.1.*x.^4-0.15.*x.^3-0.5.*x.^2-0.25.*x+1.2 


clc 

h=0.5; 

x=0.57 

x1=x+h 

fxx=[-0.1 -0.15 -0.5 -0.25 1.2] 
fx=polyval (fxx,x) 

fx1=polyval (fxx,x1) 
tr_va=polyval (polyder (fxx) ,0.5) 
fda= (fx1-fx) /h 

et= (tr_val-fda) / (tr_va1l) *100 


EXAMPLE 
Comparison of numerical derivative for backward difference and central 
difference method with true derivative and with standard deviation of 0.025 
x = [0:pi/50:pi]; 
yn = sin(x)+0.025 
True derivative=td=cos(x) 
solution 


clc 

% Comparison of numerical derivative algorithms 
x = [0:pi/50:pi]; 

length (x) ; 

% Sine signal with Gaussian random error 

yn = sin (x) +0.025*randn (1,n) ; 

% Derivative of noiseless sine signal 

td = cos (x); 

% Backward difference estimate noisy sine signal 
dynb = diff (yn) ./diff (x); 

subplot (2,1,1) 
plot(x(2:n) ,td(2:n) ,x(2:n) ,dynb, 'o') 

xlabel ('x') 

ylabel ('Derivative') 

axis([0 pi -2 2]) 

legend ('True derivative' , 'Backward difference') 
% Central difference 

dync = (yn(3:n)-yn(1:n-2))./(x(3:n)-x(1:n-2)); 
subplot (2,1,2) 
plot(x(2:n-1) ,td(2:n-1) ,x(2:n-1) ,dync,'o') 
xlabel ('x') 

ylabel ('Derivative') 

axis([0 pi -2 2]) 

legend ('True derivative', 'Central difference') 


5 
1 


Figure 1 
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Figure. Comparison of backward difference and central difference methods 


Solution 


EL 0 Feı3) 


EE Fk+1] EE و‎ 7+2] 


Example 


Consider a Divided Difference table for points following 


fae 


0.303006 


0.0000 


1.1487 


2.183 


4.9811 


(=o) TEZ tof 20,31] +(- 10)(z ج‎ 1l I, F1; 12| E 10)( BRINE t2) (o, 11,22,33] 
0.00 + ( -0.0)2.2974 + ( - 0.0)(z - 0.5)0.8418 + (x - 0.0)(z = 0.5)( - 1.0)0.3636 
2.058032 + 0.297214 + 0.36306 


disp ('*****x*** divided difference table *****x****') 
x=[]2 4 6 8 10] 
y=[]4.077 11.084 30.128 81.897 222. 62[ 
£00=y (1) 7 
for i=1:4 
£1(i)=(y(i+1)-y (1)) / (x (1+1) -x (i) ) 7 
f£O01=f£1 (1)7; 
end 
f1=[f1(1) f£1(2) f£1(3) f£1(4)] 
for i=1:3 
£2(i1)=(f£1(i+1) -£1(i))/ (x(i+2) -x(1)) 7; 
£02=f£2 (1); 
end 
f£2=[f£2(1) f£2(2) £2 (3)] 
for i=1:2 
£3(1i1)=(f£2 (i+1) -£2(i))/ (x(i+3) -x(1)) 7 
f£03=f£3 (1); 
end 
£3=[f£3(1) £3 (2)] 
disp ( IkKkkktKkktkKkkktKkktkKkkKkkKkkKkkKkKKkkKKkKkKkKkKkKKKKKKT ) 
y=input ('enter value of y:') 
p4x=£00+ ( (y-x (1) ) * £01) + ( (y-x (1) ) * (y-x (2) ) *£02+ ( (y-x (1)) * (y- 
x (2)) *£02)) 
fprintf ('\np4 (%3.3£)=%5.4f' ,y,P4x) 
syms y 
P4x=f00+ ( (y-x (1)) *f£01) + ((y-x (1) ) * (y-x (2)) *£02+ ( (y-x (1) ) * (y- 
x(2))*£02)) 


£1 = 35 9.50 25.885 70.365 
£2 = 1.5046 4.0906 11,1133 


£3 0.430 1.14 


p4x = -293/100+7007/2000*y+12037/4000* (y-2) * (y-4) 
enter value of y:8 

y3 8 

p4 (8.000)=97.3200 


Example { H.W } 
Find the divided differences (newten's Interpolating) for the data 


and compare with lagrange interpolating. 


Solution 


KkkkkkkkkkKKK divided difference table *XXXXxXxkKkkkllllk kk 
£1 = 


26. 000000000000000 4. 800000000000000 


£2 = 


-10. 600000000000000 


enter value of y:5 


p4 (5.000) =-83. 8000 
px = -283/10-53/5*y^2+419/10*y 


enter value of x:5 


Pp (5.000(=-83. 8000 
Pp = -283/10-53/5*m^2+419/10*m 


ae a ma a Sa Divided Difference table algorithm 
Sm a al a { newten's Interpolating }-------- 
clc 


disp ('******** divided difference table *********') 
x=[1 0.5 3]; 
y=]3 -10 2]; 


£00=y (1) ; 
for i=1:2; 
£1(i)=(y(i+1) -y(i)) / (x(i+1) -x(i)) ; 
f£01=f1 (1)7; 
end 
f£1=[f1(1) f£1(2)] 
for i=1; 
£2(i)=(f£1(i+1) -£1(i))/ (x(i+2) -x(i)) ; 
£02=f£2 (1) ; 
end 
£2=£2 )1( 
J5) Divided Difference table algorithm--------- 
disp (===> = { newtens Interpolating }---------—-—-——————-—-—-— 


y=input ('enter value of y:'); 

px=£00+ ( (y-x (1) ) *£01) + ( (y-x (1) ) * (y-x (2) ) * £02) ; 
fprintf ('\npx(%53.3f)=%5.4f' ,¥,PX) 

syms y 

pPx=f00+ ((y-x (1)) *£01) + ((y-x (1)) * (y-x (2)) *£02) ; 
pPx=collect (px) 


CS e mal Saa COMB aEe: WIR, sss ججج‎ 
mS ml a Lagrange's interpolation method------ 
ا ا ا ا‎ SOMBER aga دو‎ 
ا ی ا ا ا‎ Lagranges interpolation method-----------—--—-——-—- 


m=input (' enter value of x:'); 
pP=0; 

s=]1 1/2 3]; 

f=[3 -10 2]; 

n=length (s) ; 

for i=1:n; 


1=1; 
for j=1:n; 
ر (ز=~ذ) گے‎ 
1=((m-s(j))/(s(i)-s(j)))*1; 
end 
end 
pP=1.*f(i)+tp; 


end 


Pr 
fprintf ('\n p(%53.3f£)=%5.4f£' ,m,p) 


syms m 
p=0; 
for i=1:n; 
1=1; 
for j=1:n; 
£ ([=~ذ)‎ 
1=((m-s(j))/(s(i)-s(j)))*1; 
end 
end 
p=1.*f (i)+p; 
end 


p=collect (Pp) 


Example { H.W } 
Estimate the In(3) for 
Xi 


F(x) 


a) Linear Interpolation. 
B) Quardratic Interpolation 
compare between a&b 

Solution 


a) Linear Interpolation. 
F1 (x)=f (x0) + ( (Ê (x1) -Ê (x0) ) / (x1-zx0) ) * (x-x0) 
b) Quardratic Interpolation 


£2 (x) =bO+b1* (x-x0) +b2* (x-x0) * (x-x<1) 


b0= f(x0) = 0.693147180559945; 
b1= (f(x1)-f(x0))/(x1-x0) = 0.34657359027993 
b2= ((f(2)-f(x1))/(x2-x1))-b1/(x2-x0) = -0.035960259056473; 


fxl1l = 0.693147180559945-0.346573590279973 )×-2( 
inter value x:3 
fx] = 1.٢. $ A 99صش فش شڭفش9‎ 8 


fx2 = 0.346573590279973X+(-0.035960259056473X+-0.071920518112945(*)X-4( 
inter value x:3 
ÊfXx2 = 1.075681029896391 


Et1 =5.360536964281382 $% 


Et2 = 2.087293124994937 % 


Quardratic Interpolation is better than Linear Interpolation 


Vo === a) Linear Interpolation --------------------=-======-- 
o === b) Quardratic Interpolati0o n------------------------ 
Yo === compare between a&D----------------=-====-===-- 
clc 


x=input((inter value x:'); 
format long 
xi=[2 4 6]; 


disp (--------- a) Linear Interpolation----------------------- ') 
fx1=fx<(1)+((fx<2)-fx<(1))/(KiI2)-xi()))*(x-xi()) 

disp (-=------- b) Quardratic Interpolatio n----------------- ') 
bO=fx<(0); 


b1=(fx{2)-fx(1))/(i2)-xi1)); 
b2=(((fx<(§)-fx<(2))/(si(3)-xi(2)))-b1)/(xi(3)-xi(1)); 
fx2=b0+b1* (x-xi(1))+b2*(x-xi(1))*(x-xi@2)); 

% pretty(fx2)Voexpand(fx2)Vocollect(fx2) 


disp (-===------ compare between a& D------------------==-- '( 
Tv=logG); 
disp (--------- a) Linear Interpolation----------------------- ') 
Et1=abs((Tv-fx1)/Tv)*100 
disp (-===------ b) Quardratic InterpolatiOo n----------------- ') 
Et2=abs((Tv-fx2)/Tv)*100 
if Et1>E({2; 

disp('Quardratic Interpolation is better than Linear Interpolation") 
else 

disp('Linear Interpolation is better than Quardratic Interpolation') 
end 
SYMS X 
disp (-===----- a) Linear Interpolation ----------------------- 
fx1=fx<(1)+((fx<2)-fx<(1))/(KiI2)-xi()))*(x-xi()) 
disp (--------- b) Quardratic InterpolatiOo n----------------- ') 
bO=fx<(0); 


b1=(fx(2)-fx(1))/(xi(2)-xi(1)); 
b2=(((fx()-fx(2))/(xi(G)-xi(2)))-b1)/(xi(3)-xi(1)); 
fx2=bþ0+b1*(x-xi(1))+b2*(x-xi(1))*(x-Xi(2)) 


The Bisection Method for Root Approximation 


we can compute the midpoint x,, of the interval x, = x < x, with 
XxX + Kn 
= 


3 2 


ك 


Knowing x, , we can find f(x,). Then, the following decisions are made: 


1. If f(x) and f(x,) have the same sign, their product will be positive, that is, f(x) f(x) >0. 


This indicates that x, and x, are on the left side of the x-axis crossing as shown in Figure. 


In this case, we replace x, with x, . 


) ا‎ f(x,) and f(x,y, ) are 
٠ Eb n ا‎ both negative and thus 
E POET E PSE their product is positive 


. ٠ + 
KX *aN K2 1 0 X2 


Sketches to illustrate the bisection method when f(x, ) and f(x,,) have same sign 


Figure . 


2. If f(x,,) and f(x,) have opposite signs, their product will be negative, that is, f(x,,) - f(x, ) <0. 


This indicates that x, and x, are on the right side of the x-axis crossing as in Figure. In 


this case, we replace x, with x, . 


f(x, ) and f(x,,) have f(x,) and f(x,,) have 
opposite signs and thus opposite signs and thus 
their product 1s negative their product is negative 


Figure . Sketches to illustrate the bisection method when f(x, ) and f(x,,) have opposite signs 


After making the appropriate substitution, the above process is repeated until the root we are 
seeking has a specified tolerance. To terminate the iterations, we either: 


a. specify a number of iterations 


b. specify a tolerance on the error of f(x) 


Example 
Use the Bisection Method with MATLAB to approximate one of 
the roots of 


y = f(x) = 3x -2x + 6x -8 


by 

a. by specifying 16 iterations, and using a for end loop MATLAB 
program 

b. by specifying 0.00001 tolerance for f(x), and using a while end loop 
MATLAB program 


Solution: 


function y= funcbisect01 (x) ; 

FS TAN KS SDN 2S 3F OF gE =8; 

% We must not forget to type the semicolon at the end of the line 
above; 


% otherwise our script will fill the screen with values of y 


(' ددد !) طوذق 


e N ol i a aig ES EG Gi EE GE Ig a e ')( % insert line under xm and 


for k=1:16; 
f£1=funcbisect01 (x1) ; f£2=funcbisect01 (x2) ; 
xm= (x1+x2) / 2; fm=funcbisect01 (xm) ; 
fprintf ('59.6f %13.6f \n', xm, fm) % Prints xm and fm on same 
line; 

if (f1*fm<0) 

x2=xm; 
else 


grid 


Figure Î 
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لاك 


function y= funcbisect01 (x) ; 


4“ 3 ¥ 6 .F* # = 8; 


aC E. 


% We must not forget to type the semicolon at the end of the line 


above; 


% otherwise our script will fill the screen with values of y 


, xm, fm) ; 


X2=2; 


س د !) dp‏ 


while (abs (x1-x2) >2*tol) ; 
f£1=funcbisect01 (x1) ; 
f£2=funcbisect01 (x2) ; 

xm= (x1+x2) /2; 
fm=funcbisect01 (xm) ; 
fprintf ('59.6f %13.6f \n' 
if (f1*fm<0O); 


x2=xm; 

else 

x1=zxm; 

end 

end 

دد" ) 15pك‏ 
ھت دقفت % 
xm fm‏ 

1. 500000 17. 0 
1.25000 4.749023 
1.1250 1.308441 
1.06250 0.03838 
1. 0 830 -0. 506944 
1.046875 -0. 4 
1.054688 -0. 5 
1.058594 -0. 5 
1.060547 0.002604 
1.059570 -0.015168 
1.060059 -0.006289 
1. 060303 -0.0014 
1.060425 0.000380 
1. 060364 -0. 032 
1. 0604 -0.000176 
1. 060410 0.000102 


Example 


Use the Bisection Method with MATLAB to approximate one of 
the roots of (to find the roots of) 


Y=f()= x.^3-10.*x.^2+5; 


That lies in the interval ( 0.6,0.8 ) by specifying 0.00001 tolerance 
for f(x), and using a while end loop MATLAB program 


Solution: 


function y= funcbisect01 (x) ; 
y = x.^3-10.*x.^2+5; 
% We must not forget to type the semicolon at the end of the line 

above; (% otherwise our script will fill the screen with values of y) 


while (abs (x1-x2) >2*tol) ; 
f£1=funcbisect01 (x1) ; 
f£2=funcbisect01 (x2) ; 
xm= (x1+x2) /2; 

fm=funcbisect01 (xm) ; 
fprintf ('%59.6f %13.6f \n', xm,fm) 7; 
if (f1*fm<0O); 


x2=xm; 

else 

x1=xm 

end 

end 

۴F)‏ ا چ ي چ ي اي ا ا س ا ا ا 
کک تف قفتت تف تققد دد % 
xm fm‏ 
0.440 0.700000 
5-- 0.750000 
0.12488 0.725000 
2 - -- 0.737500 
0.04753 0.731250 
7 :0.734375 
53- 0.735938 
0.007228- 0.735156 
0- 0.734766 
0.004 0.734570 
0.000843- 0.734668 
0.04- 0.734619 
5 :0.734595 
0.05- 0.734607 


Newton-Raphson Method 


The Newton-Raphson formula can be derived from the Taylor series expansion 
of f(x) about x: 
fn) = fOD + f D(x — x?) + Ox — XP (a) 
Ifx;,ı is a root of f(x) = 0, Eq. (a) becomes 
0= f(x) + fF (xD (x1 — xı) + O(Xr1 — xX) (b) 
Assuming that x; is a close to X;,ر,‎ we can drop the last term in Eq. (b) and solve for 
X;,1. The result is the Newton-Raphson formula 


f) 
FO) 
Ifx denotes the true value ofthe root, the error in x; is E; = xX — X;. Itcan be shown 
that if x, is computed from Eq. ( c ), the corresponding error is 
(xD ~2 
2F 


ams" 
le 


Xi] = Xj — 


E1 = 


indicating that the Newton-Raphson method converges quadratically (the error is the 
square of the error in the previous step). As a consequence, the number of significant 
figures is roughly doubled in every iteration, provided that x; is close to the root. 


f(x) 
Figure a )Graphical interpretation of the Newton-Raphson 
formula. 


X1 


A graphical depiction of the Newton-Raphson formula is shown in Fig. ( a ) The for- 
mula approximates f(x) by the straight line that is tangent to the curve at X;. ThUS Xj, 1 
is at the intersection of the x-axis and the tangent line. 

The algorithm for the Newton-Raphson method is simple: it repeatedly applies 
Eq. (c ), starting with an initial value xo, until the convergence criterion 


xı — xı| < € 


is reached, £ being the error tolerance. Only the latest value of x has to be stored. Here 
is the algorithm: 


1. Let x be a guess for the root of f(x) = 0. 
2. Compute Ax = - f(X/ F(X). 
3. Letx > x + Ax and repeat steps 2-3 until |AX| < #. 


EXAMPLE 


A root of f(x) = x3 - 10x + 5 = 0 ¡lies close to x = 0.7. Compute this root with the 
Newton-Raphson method. 


Solution 


The derivative of the function is f(x) = 3x7 —- 20x, so that the Newton- 
Raphson formula in Eq. ( c ) is 
E REE x32 1025 21 10225 
f 3x2 - 20x, x%(3x- 20( 
It takes only two iterations to reach five decimal place accuracy: 
20.7) - 10)0.7(2 5 a 
Pp i eee CL? Ff 
“ 0.780.7 - 20| 
+ ے‎ 5 36) - 10)0.735 36) - 5 
0.735 36 ]3)0.735 36) - 20[ 


= 0.734 60 


Example 


Use the Newton-Raphson Method to estimate the root of f(x<)=e^(-x)-x, 
employing an initial guess of x0=0 


x=[0]; 
tol=0. 0000000007; 
format long 
for i=1:5; 
fx=exp (-x(i1i)) -x(3i) ; 
fxx=-exp (-x(i1)) -1 ; 
fxxx=exp (-x (1) ) ; 
x (1+1) =x (i) - (fx/Êxx) ; 
T.V(i)=(abs ((x(i+1) -x(i) ) /x (1+1) ) ) *100; 


end 
for i=1:5; 
e(i)=x (6) -x(i) 7 
fxx=-exp (-x(6))-1 ; 
fxxx=exp (-x (6) ) ; 
e(i+1)=(-fxxx/2*fxx) * (e (i) ) ^2; 
end 


if abs (x(i+1) -x(i) ) <tol 
disp (' enough to here') 
15) '( 
disp(' X(i+1) ') 
ISP (SSS SSE '( 


enough to here 


0.500000000000000 
0.566311003197218 
0.567143165034862 
0.567143290409781 
0.567143290409784 


1.0e+002 * 


1.000000000000000 
0.117092909766624 
0.00114 5 
0.0000002210619 
0.000000000000005 


0.567143290409784 
0.067143290409784 
0.000827256 
0.000000125374922 
0.000000000000003 
0.000000000000000 


The secant Formula Method 


A popular method of hand computation is the secant formula where the improved 
estimate of the root (x,,,) is obtained by linear interpolation based two previous 
estimates (xX, and Xx-1): 


JA = XI -1 


e a A RE 3 1 
fx) — f(1) 


J1 = Xj = 


Example 


Use the The secant Formula Method to estimate the root of 
f(x)=e^(-x)-x, employing an initial guess of x(1-1)=0 & x(0)=0 


A = X-1 


I] 5 J = 1 (x; 1 
FD = 
Solution 
ا ڪاڪ ڪڪ .۽‎ The secant Formula Method ==------------= === 
clc 
x=[0 1]; 


TV=0.567143290409784; 
format long 
for i=2:6; 
fx=exp (-x (i-1) ) -x(i-1) ; 
fxx=exp (-x(i)) -x(1) ; 
x(i+1)=x (i) - ((x (i) -x (i-1) ) *fÊxx) / (fxx-fÊx) ; 
E_T(i)=(abs ( (TV-x(i+1) ) /TV) ) *100; 
end 


. 000000000000000 
. 612699836782 
. =¬ 4 
. 56717035845 
. 567143306604963 
. 5671432904975 


OOOOO 


N 8 
. -. 2 
N 40 
. 00000285557096 
. 000000000017 


OOOO 0O 


Example 


Use N.R. Quadratically Method to estimate the multiple root of 
f(x)=x^3-5x^2+7x-3, Initial guess of x(0)=0 


f) F(0) 


Xr = ل‎ = 
+1 1 2 


f(u— f(x) FC 


Solution 
=== The N.R. Quadratically Method  ------—-———————-—-—- 
clc 
TV=1; 
x= [0]; 
format long 
for i=1:6; 
fx=x (1) ^3-5*x (i) ^2+7*x (1) -3 
fxx=3*x (1) ^2-10*x (1) +7 
fxxx=6*x (1) -10 
x (1+1) =x (i) - (fx*fÊxx) / ( (fxx) ^2-fx*fÊxxx) 7; 
E_T(i)=(abs ( (TV-x(i+1) ) /TV) ) *100; 
end 
dISp (SSeS ۴ 
disp(' X(i+1) ') 
313p) جد ددد‎ '( 
× ۲ 
IIS (aaa '( 
disp (' ET ')( 
dISp( Seas '( 
E_T' 
disp (== '( 
ك‎ e e 
ڪڪ‎ Multiple Roots-------—--—- 
$--fx= (x-3) (x-1)( (x-1( --------- 
clc 


for x=-1:0.01:6; 
fx=x.^3-5.*x.^2+7.*x-3 
plot (x, fx) 
hold on 

end 

grid 

title (' (x-3) (x-1) (x-1) ') 

xlabel ('x') 

ylabel (' fx') 


1.113635789 7 
1.003081664098603 
1.000002384 6 
1.00000000007312 
1.00000000007425 
1.000000000074625 


1.3363436 
0.308166403 
0.00023814938 
0.000005 
0.00000000746275 
0.000000075 


Figure Î 


File Edt View Inset Tools Desktop Window Help 
TLE IEEWE OEE 
[ Fx=(-Sx-x-1) ] 


20 


BU 
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Example 


Use the Newton-Raphson Method to estimate the root of 
f(x)=x^3-5x^2+7x-3, Initial guess of x(0)=4 


x=[4]; 

tol=0.0007; 

TV=3; 

format long 

for i=1:5; 
fx=x (i1) ^3-5*x (i) ^2+7*x (1) -3; 
fxx=3*x (1) ^2-10*x (i) +7; 
x (1+1) =x (i) - (fx/Êxx) 7 
E_T(i)=(abs ( (TV-x (i+1) ) /TV) ) *100; 

end 

for i=1:5; 


e(i)=x (6) -x(i) 7 
fx=x (1) ^3-5*x (1) ^2+7*x (1) -3; 
fxx=3*x (1) ^2-10*x (i) +7; 
fxxx=6*x (1) -10; 
e(i+1)=(-fxxx/2*fÊxx) * (e (i) ) ^2; 
end 
if abs (TV-x (i+1) ) <tol 
disp(' enough to here') 


enough to here 


4. 000000000000000 
3. 400000000000000 
3. 100000000000000 
3.00866213 
3. 00007492 
3. 000000005570623 


CEKCKEKEEEEKEKEEEEEEI] 
CEKCKEKEEEKEEEEEEYY 
0.289855072463781 
0.0024880263730 
0.000000185687436 
0.000000007462475 


-0. 999999994427 
-0.39999999447 
-0.099999994427 
-0.00869564660 
-0.0000746352 569 
-0.00000008914454 


ترقبوا المزيد من الشروحات للأمثلة فى التحليل 
العددى والرياضيات والتحكم الألى والأتصالات 
وإلكترونات القدرة ونظم التشغيل والدارات التماثلية 
والنظم الرقمية واسس الألكترونات وغيرهامن 
المواد فى اغلب التخصصات راجين من الله سبحانه 
وتعالى التوفيق فلله الحمد والمنة وبه التوفيق 
والعصمة. 


